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ANALYTICAL METHOD FOR DETERMINING THE STABILITY 
OF LINEAR RETARDED SYSTEMS WITH TWO DELAYS 

L. Keith Barker 
Langley Research Center 

J. L. Whitesides 

The George Washington University 
Joint Institute for Acoustics and Flight Sciences 

SUMMARY 

Differential -difference equations occur in various branches of science. The sta- 
bility of these equations, however, is much more difficult to determine than that of ordi- 
nary differential equations. In this paper the stability of differential-difference equations 
of the retarded t 3 q)e with constant coefficients and two constant time delays is considered. 
A new method that makes use of analytical expressions to determine stability boundaries, 
and hence the stability of the equations, is derived. The basis of the method consists in 
deriving analytical equations for each of the delays which correspond to the purely imag- 
inary roots of the characteristic quasi -polynomial. 

The method developed is used to analyze the stability of a second -order differential 
equation with delays in the velocity and displacement terms. The resulting stability re- 
gions are in agreement with those obtained by other investigators. 

INTRODUCTION 

A differential -difference equation of the retarded type is basically a differential 
equation in which the highest order derivative of the dependent variable contains no delay 
in its argument (time), but any of its other derivatives or the dependent variable itself 
may have delays. Such equations arise when the future state of a system depends not 
only on its present state but also on part of its past history. (See refs. 1 to 5.) 

One of the major obstacles encountered in dealing with systems with delays is the 
stability analysis. Moreover, it appears that practical applications of such systems have 
been limited thus far to one and two delays. A convenient method is developed in refer- 
ence 6 to examine the stability in the delay space for a linear, time invariant differential - 
difference equation of the retarded type with one constant delay. This method, which is 
somewhat similar to the method of reference 7, is extended in reference 8 to many con- 


stant delays. In particular, stability boundaries in the delay space for a retarded system 
with two delays can be constructed by using the method of reference 8. However, the 
points on the boundaries are found by an indirect procedure which involves searching for 
the roots of transcendental equations, (This is not the case for one delay.) 

The purpose of the present study is to develop explicit analytical equations for con- 
structing stability boundaries in the delay space for retarded systems with two constant 
delays. This report is based on work that has been completed as part of L. Keith 
Barker’s doctoral research program with The George Washington University. 

SYMBOLS 

A, A; N X N matrices of real constants 

a,b,c real constants 

^lh^l2>^21>^22 constants 

B N X m matrix of real constants 

constants 

C positive scalar function of co (see eq. (43)) 

I N X N identity matrix 

i imaginary unit, 

i,k,Z,m,n integers 
Kj integer associated with 

L(s),L(s,0j) characteristic quasi -polynomial of retarded system 
Lq(s) resulting polynomial when delays are made zero in L(s) 

M number of delays in system 
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N dimension of system 

P(s),Q(s),R(s) polynomials in s with real constant coefficients 
q,r angles defined in figure 1 

s complex variable, a + io) 

s* root of L^s,0jj 

t time or independent variable 

t nondimensional time 

if m X 1 vector forcing function 

Vj integer associated with 

X scalar function of time 

X N X 1 state vector 

yi,y 2 defined by equations (14) and (15), respectively 

a argument of Q(ict)) (see fig. (1)) 

/3 argument of P(iw) (see fig. (1)) 

y argument of R(io)) e ^^^2 (ggg (j)) 

5,6 small positive numbers 

C damping parameter 

real and constant delays 

p radius of small circular contour around any root of L(s,0j) 
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a 


real part of complex variable s 



real and constant time delays 


imaginary part of complex variable s 


largest positive real root of equation (20) 

Mathematical notations : 

11 

absolute values or magnitude 

arg 

argument 

det 

determinant 

min 

minimum value 

A 

small incremental value 


Dots over a symbol denote derivatives with respect to time. 

ANALYSIS 
Retarded Systems 

The state equation for a physical system is sometimes written in the form 

M 

x(t) = A x(t) + ^ A^ X (t - T^) + B u(t) (1) 

f=l 

where x (t) is an N x 1 state vector, u(t) is an m x 1 input vector, A and A^ 
are N x N constant matrices, B is an N x m constant matrix, and = 0 is a 

constant time delay. Equation (1) is a linear differential-difference equation of the re- 
tarded type with constant coefficients and delays. For brevity, equation (1) is referred 
to as a retarded system. Mathematical treatments of retarded systems can be found in 
references 1, 9, and 10, for example. 
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Equation (1) is said to be stable if and only if the corresponding homogeneous 
equation 


M 

X (t) = A x(t) + ^ X (t - r-j^ 
l=\ 


(2) 


is stable; that is, if all the roots s = a + io) of the characteristic quasi -polynomial 


L(s) = det 


si 




( 3 ) 


have negative real parts. In this case, the solution to the homogeneous equation (2) 
approaches zero as t — ■» (ref. 1). 

Techniques of generalized harmonic analyses can be applied meaningfully to the 
retarded system of equation (1) when equation (2) is stable (ref. 11). The stability of 
equation (2) is also inherently related to the asymptotic stability of some nonlinear con- 
trol problems (ref. 10). Various methods have been proposed for determining whether 
any of the roots have nonnegative real parts (cr = 0); for example, see references 8, 12, 
13, and 14. 


Class of Retarded Systems 

Stability boundaries are examined for the class of retarded systems which has a 
characteristic quasi -polynomial of the form: 


L(s) = P(s) + Q(s) e'^^^ + R(s) e (4) 

where 

N 

P(s) = ^ ajs^ (ajvf # 0) 

j=0 


Q(s) 


N-1 
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R(s)= ^ Cjsj 
j=0 

The real nonnegative constants and 62 will be referred to as delays although 
they may be some linear combination of the actual time delays in equation (1), Some 
examples of retarded systems with characteristic quasi -polynomials of the form of equa- 
tion (4) are given in appendix A. 

The stability condition (stable or unstable) of a retarded system with characteristic 
quasi -polynomial equation (4) is the same as the same system with zero delays 
(0j = 02 = 0), if the values of Q-^ and 02 equation (4) are sufficiently small. This 
initial stability condition is determined by setting 0^ = 02 = 0 in equation (4) and then 
examining the roots of the resulting polynomial 

Lo(s) = P(s) + Q(s) + R(s) (5) 

If the degree of the polynomial equation (5) is N and if the delays S\ and 02 are 
sufficiently small, then the quasi -polynomial equation (4) has N roots very close to the 
N roots of equation (5), and the remaining infinity of roots have very large negative real 
parts (refs. 8 and 13). Additional information on the roots of quasi -polynomials is con- 
tained in reference 4 and 15. 

As the delays are varied in some continuous manner from essentially zero, the 
roots of the quasi -polynomial equation (4) move continuously and generate an infinite 
number of continuous root -locus curves in the complex root plane; that is, in the cto>- 
plane, s = ct + iw satisfies L(s) = 0. (See appendix B.) Clearly, it is impossible to 
plot all the root -locus curves of equation (4) for a given continuous variation of the de- 
lays. It is possible, however, to determine the number of roots with positive real parts 
as the delays are varied by examining the behavior of the root -locus curves on the imag- 
inary axis. This approach is used in reference 8. 

Partitioning Delay Space Into Stable or Unstable Regions 

Root-locus curves are generated by the roots of equation (4) as the delays are 
varied. If a root -locus curve comes into contact with, or crosses, the imaginary axis 
at w (touch point), then by definition, s = icu satisfies 

L(io)) = P(iw) + Q(ia>) e + R(io;) = 0 (6) 
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Since the root -locus curves are symmetrical about the real axis, only w ^ 0 values are 
considered. 

Equation (6) can be considered as two equations (real and imaginary parts) in two 
unknowns. The two unknowns which are chosen vary with the different methods which 
have been used to examine the stability of systems with delays. For example, in refer- 
ence 6, 02 = 0 and the two unknowns are w and in Neimark's method (ref. 14), 
the delays are held fixed, and the two coefficients are chosen as the unknowns with <x> 
as a coordinating parameter; in reference 16, the unknowns are a delay and a coefficient 
or gain. In the present study, the two unknowns are d\ and 02 with w as a coor- 
dinating parameter. 

A vector representation of the complex quantities appearing in equation (6) is shown 
in figure 1. The two distinct solution sets to equation (6) are shown graphically in fig- 
ure 2. From the geometric properties of a triangle, it follows that a solution to equa- 
tion (6) exists if and only if the following three relationships simultaneously hold: 


|P(ico)| + lQ(io))| £ |R(iw)| 

(7) 

|P(M| + |R(iO))| ^ |Q(ico)( 

(8) 

|R(iw)| + Q(iaj) s |p(io))| 

(9) 


These relations express the fact that the sum of the lengths of any two sides of a triangle 
must be greater than or equal to the length of the remaining side. An equality sign in 
either of equations (7), (8), or (9) corresponds to collinear vectors. 

It follows by using figure 1 that the angles r and q in figure 2(a) are given by 
r = 77 - /3 + 0 ! 

r = 77 - arg P(ico) + arg Q(ico) - u)9\ (10) 


and 


q = 77 - ^ + y 

q = 77 - arg P(ici>) + arg R(iw) - co 6 2 (11) 
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Now, solving equation (10) for 6 i and equation (11) for 62 gives 


9i = -^[77 - arg P(io)) + arg Q(iw) - r] 


(: 


and 


02 = - arg P(io)) + arg R(ia>) - 


The angles r and q are obtained from figure 2(a) by applying the law of cosines as 


cos r 


_ |P(M|^+ |Q(la))|2 - |R(ia))|2 _ 

2|P(ico)||Q(ia))| 


(] 


and 


|p(M| 2 + |R(ico)|2 . |Q(io,)|2 

COS (-q) = ^ ^ ^ = Vo 

2|P(iw)||R(io))| ^ 


Choosing 0 ^ cos"^ yj ^ tt and using the geometry in figure 2(a) gives 
Itt - arg P(ia>) + arg Q(iw) - cos"^ y^ + 2irK^ 


01 = i 

CO 


and 

7T - arg P(ico) + arg R(ico) + cos'^ y 2 + 277 X 2 ' 

where Ki and K 2 are integers. 

Using the geometry in figure 2(b) gives 



01 = 1 77 - arg P(ico) + arg Q(ico) + cos'^ y, + 2 t7Vi1 


02 = - arg P(io)) + arg R(ico) - cos"^ Y 2 + 27rV2j 


(19) 


where Vj and V 2 are integers. A solution to equation (6) exists if and only if 0j 
and 02 satisfy the pair of equations (16) and (17) or the pair of equations (18) and (19). 

Combination values of 0j and 02 for which a root -locus curve touches the 
imaginary axis can be determined by plotting 0j and 02 against co. Alternately, 02 
can be plotted directly against 0 1 for corresponding values of w. This latter type of 
figure is actually a partitioning of the delay space into regions of stability (stable or un- 
stable). It should be noted that Neimark's D-partition methods (ref. 14) can be used to 
construct regions of stability in the plane of two real parameters (gains or coefficients) 
which occur linearly in the characteristic quasi -polynomial for fixed delays. Other 
methods, developed within the last decade, are discussed briefly in reference 17. 

The stability region for the smallest values of delays may be all that is required 
to show that the delays are not large enough to make the system unstable or to indicate 
how much the delays can be increased before the system becomes unstable. 

Special Values of a> 

There are certain values of co which are useful in evaluating the equations for 
0j and 02 > 

Upper bound on o) . - An upper bound on o) can be computed by using equation (6). 
The dominant power on w occurs, by definition, in P(io)); hence, let w = com be the 
largest positive real root of the equation 



N-1 


^ (l^nl 


+ Ibr 


|Cn|)o)« = 


= 0 


n=0 


( 20 ) 


Then, it follows that co S com in equation (6) and in the pertinent pair of equations for 
01 and 02. 

Border values of ct) .- Partitioning curves are defined only for those values of o) 
which satisfy equations (7) to (9). These meaningful values of o) are determined by 
using border values of o). A border value of w is defined as a nonnegative real value 
of w which satisfies any of the equality relations in equations (7) to (9), which are: 


|P(iw)| + |Q(io))| = |R(io))| 


(21) 
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lP(iw)l + |R(iw)l = |Q(io))l 
|R(io))| + |Q(io.)| = |P(iw)| 


( 22 ) 


(23) 

These equality relationships also follow by setting cos r = ±1 or cos (-q) = ±1 in equa- 
tion (14) or (15), respectively. 

The finite number of border values separates the w space into different intervals. 
To determine if a partitioning curve is defined for values of w in an interval, evaluate 
cos r or cos (-q) at some value of w within the interval. If, for example, |cos r| 2 1, 
then partitioning curves exist for all values of w in the interval. On the other hand, if 
I cos r| > 1, the partitioning curves do not exist for any value of w in the interval. The 
maximum border value will be equal to cum. 


Multiple Touch Points 

Multiple touch points occur when more than one root -locus curve touches the 
imaginary axis for the same values of oj, 6i, and 02« this case, the points on the 
partitioning curves correspond to purely imaginary roots of the characteristic quasi - 
polynomial which have multiplicities greater than imity. Because of its intuitive appeal, 
a procedure used in reference 8 is adapted later to examine the stability condition of 
regions established by the partitioning curves. This procedure, however, is presented 
in reference 8 and in this paper on the assumption that the touch point, where the proce- 
dure is applied, is simple. This condition is easily checked by using the partitioning 
equations. 

By construction, s = icu is a root of the characteristic quasi -polynomial for any 
point (01,02) the delay space which lies on a partitioning curve. The questions is: Is 
s = io) a simple or multiple root? 

Differentiating equation (4) with respect to s gives 


ds ds 







(24) 


The value of dL/ds at a point (0i,02) on a partitioning curve is given by equation (24) 
with s = i(^ and with 0i and 02 computed by using the partitioning equations. If 
dL/ds it 0 at this point, then s = i(^ is a simple root; otherwise, it is not. 

Stable and Unstable Regions in Delay Space 

The delay space is partitioned into regions which are either stable or unstable by 
the 02 versus 0i. The stability of each of these regions can be determined by com- 
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puting the initial stability of the system using equation (5) and evaluating the derivative, 
on a partitioning curve, of the real part of s = a + ico with respect to one of the delays 
while holding the other delay constant. 

Select one of the delays 6i or 02 to be held fixed. For discussion purposes, 
denote this delay by 0; (i = 1 or 2) and the remaining delay by 0^ (k i; k = 1 or 2). 

It is assumed that the purely imaginary root is simple (dL/ds o) at a point on the 
partitioning curves. Moreover, suppose that at this point a^a/a0ij^ = 0 for 
n = 0, 1, . . ., g - 1; but that d^a/d9\^^ 0. Then, if d^a/d9y^ > 0 and g is an odd 

integer, the sytem gains a root with positive real part when the point (0k.^i) is varied 
from a partitioning curve in the direction of increasing 0jj. Likewise, if d^(j/d9^ < 0 
and g is an odd integer, the system loses a root with positive real part. If 
a^a/a0jjS ii o and g is an even integer, then the root-locus curves are tangent to the 
imaginary axis, so that the system neither gains nor loses any roots with positive real 
parts. 

For convenience the characteristic quasi -polynomial is written as 

L = L(s,0^,0k) (25) 

The characteristic equation is then 

L = L(s,0z,0k) = 0 (26) 


Applying implicit differentiation to equation (26) and extracting the real part of as/a0ij 
gives 


where the partial derivatives = aL/90k and Ls = aL/9s are obtained by using 
equation (25) and R denotes the real part of the complex number in parentheses. 

The derivative in equation (27) is evaluated on a partitioning curve by setting 
s = ico in equation (27) and using the partitioning equations to obtain 9i and 0 k. This 
derivative then gives the change in the real part of the characteristic root when leaving a 
partitioning curve at (0i,0k) iii the direction of increasing 0k. 

If 9a/90k = 0, then a similar procedure is followed for higher derivatives until 

d^ajd9^ 0 . 
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Summary of Stability Procedure 

The procedure used to partition the delay space into different regions and to iden- 
tify the different regions as stable or unstable is summarized in the following steps: 

(1) Compute the initial stability of the system by using equation (5). 

(2) Compute meaningful values of co by using equations (21) to (23). 

(3) Plot partitioning curves in the delay space by using 6i and 02 equations as 
w varies over the predetermined range of acceptable values in step 2. 

(4) If the delays do not cross a partitioning curve (in the plot of step 3) as they are 
varied in combination from zero to their final desired values, then the retarded system 
with characteristic quasi -polynomial of the form of equation (4) maintains its initial sta- 
bility. However, if the delays cannot be varied to their final values without crossing a 
partitioning curve, compute the multiplicity of the touch points and the derivative of the 
real part of s = a + iw as the curve is crossed. Use equations (24) and (27). 

The stability on each side of a partitioning curve is determined by counting the 
number of roots with positive real parts as the curves are crossed. Stability boundaries 
are those curves which partition the delay space into stable and unstable regions. 

EXAMPLE 

To clarify the stability analysis presented in the previous section, the following 
example is examined: 

i^x(t) + 2?4x(t -0i) +x(t - 02 ) = O (28) 

dt^ dt 

where ? > 0 is a real constant damping parameter, and the delays 0^ and 02 are 
nondimensional. This example is examined in reference 8 in a different manner. 

The characteristic quasi -polynomial for equation (28) is 


L(s) = P(s) + Q(s) e + R(s) e”^2S 


(29) 


where 


P(s) = s2 Q(s) = 2?s R(s) = 1 
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The initial stability of the retarded system is found by using equation (5), which is 


Lq(s) = s^ + 2?s + 1 


(30) 


Since ? > 0, the roots of Lq(s) have negative real parts, so that the retarded system 
is initially stable. 

Setting s = io) in equation (29) gives 
P(iw) = 

Q(iu)) = 2?0)i 


R(iw) = 1 


The first pair of equations for the delays Q\ and 02 (eqs. (16) and (17)) becomes 


01=1 


7T 

^ - COS 
2 


-1 


4Co)' 


^ + 2vKi 


(31) 


and 



cos 


-1 



- 


2w^ 



+ 2itK2 


(32) 


The second pair of equations for 0i and 02 (eqs. (18) and (19)) becomes 


H=h 


^ + cos 

a 


-1 - i\ 


4?cd‘ 


27tVi 


(33) 


and 


0o = 1 
<5 O) 


_1 (( j )^ - 4?2co2 +1^ 

-cos ^ ( — 2 j + 2V V2 


(34) 
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An Upper Bound and Border Values of co 


An upper bound Wm.- The upper bound of o) is obtained by using equation (20), 
which becomes 

0)2 - - 1 = 0 (35) 

The largest real nonnegative value of o) which satisfies equation (35) is denoted by o)m. 
In this case, 

= ? + \ZC^ + 1 (36) 


Border values of o) . - The border values of o) are obtained by using the relations 
in equations (21), (22), and (23), which become 


0)2 + 2 ^ 0 ) -1 = 0 


(37) 


0)^ - 2?o) +1 = 0 


(38) 


- 2 ^ 0 ) -1 = 0 


(39) 


Notice that equation (39) is the same as equation (35), so that 0)^^ is a border value. 
Other border values are obtained from the solutions of equations (37) and (38) which are, 
respectively, 



(40) 

= ?±\/?2 _ 1 

(41) 


Only the positive radical in equation (40) is of interest since o) 5 0. Equation (41) is 
only of interest when ? = 1. 


Specific Calculations 

The border values and upper bound of co are shown in table I for damping param- 
eter values ? of 0.2, 0.5, and 1. 
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TABLE I.- BORDER VALUES AND UPPER BOUND OF w 


Damping parameter, ^ 

Border values of ct) 

Upper bound, ctim 

0.2 

0.82 

1.22 


1.22 


.5 

.62 

1.62 


1.62 


1.0 

.41 

2.41 


1.00 



2.41 



Figures 3 and 4 show the pairs of the delays and 02 which result in a root- 
locus curve coming into contact with, or crossing, the imaginary axis in the complex root 
plane at co (touch point) for various values of (Ki,K 2 ) and (Vi,V 2 ), respectively. Re- 
sults are presented for ^ = 0.2, 0.5, and 1. 

The consecutive Kj curves in figure 3 and the consecutive Vj curves in fig- 
ure 4 both differ by 2jr/o>, as shown by equations (31) to (34). Thus, the V 2 = 0 curve 
falls off the scale in figure 4. 

The terminal points of the curves in figures 3 and 4 correspond to a zero delay or 
to a border value in table I. All border values, however, do not necessarily specify a 
terminal point on the curves. The value w = 1 is a border value that occurs along the 
curves in figures 3(c) and 4(c) and is identified by a singularity in the slope of the curves 
at this value. 

Figures 3 and 4 can be used to obtain the touch points cu which occur as the delays 
are varied in some continuous manner from zero to their final constant values. The 
touch points clearly depend on the manner in which the delays are varied. 

A partitioning of the delay space results when 02 is plotted directly against 0j 
with CO as a coordinating parameter, as shown in figure 5. The solid curves correspond 
to various values of (Ki,K 2 ) and are generated by using the pair of equations (31) and (32). 
The dashed curves correspond to various values of (Vi,V2) and are generated by using the 
pair of equations (33) and (34). The totality of curves for (Ki,K 2 ) and (Vi,V 2 ) partition 
the delay space into regions. The arrows on the curves denote the direction of increas- 
ing Cl). 

Since the system is initially stable, it will remain stable until the curve 
(Ki,K 2) = (0,0) In figure 5 is touched. To examine the stability beyond this point re- 
quires the determination of the multiplicities of the touch points and the changes in the 
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real parts of the root -locus curves as they pass beyond the touch point. All the regions 
shown in figure 5 can be entered by crossing one of the solid curves corresponding to 
(Ki,K 2); hence, only these curves need to be examined further. 

It can be shown by setting s = ioj in equation (24) and using equations (29), (31), 
and (32) that dL/ds 0 on any of the (Ki,K2) solid curves of figure 5, Hence, the mul- 
tiplicity of all touch points is unity. 

The change in the stability at a touch point is examined by using equation (27) with 
0^ = and Qi = 02- Reference 8 expresses equation (27) for this particular application 
in a very simple form which is equivalent to 


= 0 )“^ + 0)202 0 ) 02 - 1 

90 


(42) 


where 


C 


2 ^ 


1 

9L(iCL)) 


901 


(43) 


The scaled derivative (90/90 j) has the same sign as 9a/90i. The scaled derivative 
on the (Ki,K 2) curves in figure 5 is evaluated by using equation (32) to obtain 02. Since 
9 1 does not appear in equation (42), the scaled derivative is the same for all curves 
(Kx,K 2) which have the same K2 value. The scaled derivative is plotted in figure 6 as 
a function of 02 for three values of K2, which are 0, 1, and 2. The arrows on the 
curves indicate the direction of increasing o). The K2 = 0 curve in figure 6 applies to 
the (0,0), (1,0), (2,0), and (3,0) curves in figure 5; the K2 = 1 curve applies to the (0,1), 
(1,1), and (2,1) curves; and the K2 = 2 curve applies to the (0,2), (1,2), and (2,2) curves. 

Figure 6 is used to determine the stability condition on each side of the partitioning 
curves in figure 5. The arrows on the curves denote increasing values of w and are 
used for corresponding values of the scaled derivative with points on the partitioning 
curves. For instance, suppose 02 = 0,2 on the (Kj,K 2) = (0,0) curve in figure 5, then 
the corresponding value of scaled derivative in figure 6 is 0.98. This correspondence is 
established by noting that 02 decreases from 02 = 0.2 with increasing o). Now, 
since C^[da/d9i)>0 on (Ki,K 2) = (0,0) when 02 = 0.2, it follows that proceeding off 
the (Kx,K 2) = (0,0) curve in figure 5, with 02 = 0.2, in the direction of 0J, results in 
the gain of a root with positive real part. The system, therefore, becomes unstable. 

Next, consider the curve (Kj,K2) = (1,0) in figure 5. It can be shown by using 
figure 6 that entering the lower triangular region from the left results in the loss of a 
root with positive real part, so that the system becomes stable again. Upon leaving the 
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triangular region on the right, however, a root with positive real part is gained, and the 
system becomes unstable. 

Entering the upper looped portion of the curve (Ki,K 2 ) = (1,0) from the left 
results in the gain of a root with positive real part. So, inside the upper looped region, 
the system has two roots with positive real parts. Leaving the upper looped portion of 
the curve on the right results in the loss of a root with positive real part; however, the 
system had two roots with positive real parts inside the looped portion of the curve, and 
is, therefore, still unstable. 

This type of reasoning is repeated for the other regions of figure 5. It is found 
that there are only two stable regions in figure 5. These are the initial stability 
region bounded by the curve (Ki,K 2 ^ = (0,0) and the lower triangular region of the 
(Ki,K 2 ) = (1,0) curve. 

Figure 7 shows the stability boundaries for the retarded system represented by 
equation (28) when ^ = 0,2, 0,5, and 1. The hatched lines indicate the stable side of the 
boundaries. These results are in agreement with those of reference 8. 

It is well known that the solution of equation (28) becomes more highly damped as 
^ increases if 61 = 62 = 0 . If 6 i* 0, this may not be the case, as shown in figure 7. 
For example, let 62 = 0 and 6 i = 0.8. Then, the solution of equation (28) is stable for 
^ = 0.2 and unstable for ? = 1. 

CONCLUDING REMARKS 

A new method is developed for generating stability boundaries for differential- 
difference equations of the retarded type with constant coefficients and two constant de- 
lays. The basis of the method consists in deriving analytical equations for each of the 
delays which correspond to the purely imaginary roots of the characteristic quasi - 
polynomial. These analytical equations are then used to partition the delay space into 
regions which are stable or unstable. 

Each point on a partitioning curve corresponds to a root -locus curve touching or 
crossing the imaginary axis in the complex root plane. The stability of the different re- 
gions, therefore, is examined by computing the initial stability of the equations (zero 
delays) and counting the number of root -locus curves with positive real parts as the dif- 
ferent regions are entered. The analytical equations mentioned previously are used in 
counting the number of root-locus curves with positive real parts. Specifically, the ana- 
l 5 d:ical equations are used in computing the multiplicities of the root -locus curves and the 
partial derivative of the real part of the root -locus curves with respect to one of the 
delays. 
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The method developed is used to analyze the stability of a second -order differential 
equation with delays in the velocity and displacement terms. The resulting stability 
regions are in agreement with those obtained by other investigators. 

Langley Research Center 

National Aeronautics and Space Administration 
Hampton, Va. 23665 
September 12, 1975 
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APPENDIX A 


EXAMPLES OF RETARDED SYSTEMS WITH TWO DELAYS 


The followir^ examples show retarded systems which have characteristic quasi - 
polynomial equations of the form of equation (4): 

Example 1: Consider the following scalar differential equation with two constant real 
time delays tj g 0 and T 2 = 0: 

N N-1 N-1 

^ anX^N)(t) + ^ bnx^”\t - t^) + ^ c„ x^*^\t - xg) = f(t) (Al) 

n=0 n=0 n=0 

where x^*^^(t) denotes the nth derivative of x(t) and f(t) is a continuous input func- 
tion of time. The coefficients an, bn, and Cn are real constants. Equation (Al) can 
be written in the form of equation (1); however, this is unnecessary in obtaining the 
characteristic quasi -polynomial. The characteristic quasi -polynomial associated with 
equation (Al) is equation (4) with 


N 

P(s) = ^ ans" 
n=0 


Q(s) 


N-1 


n=0 


bns« 


N-1 
R(s) = 2 
n=0 


CnS‘ 


01 = T1 


02 = T2 
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Example 2: Let the retarded system be described by two coupled linear equations with 
one distinct delay as follows: 


Xl(t) 


^11 ^12 

Xi(t) 

+ 

bll bi2 

xi(t - t) 

+ 

fl(t) 

X2(t) 


^21 ^22_ 

X2(t) 


•^21 b22_ 

X2(t - t) 


^2(^) 


The characteristic quasi -polynomial associated with equation (A2) is given again by 
equation (4) with 

P(s) = s2 - (ail + ^22)® + ^11^22 ‘ ^21^12 

Q(s) = -(bii + b22)s + aiib22 + a22t>ii - a2ibi2 - ai2b2i 

R(s) = biib22 - b2lbi2 

01 = T 

02 = 2r 


Note that although there is only one time delay t in equation (A2), there are two delays 
01 and 02 in the associated characteristic quasi -polynomial. 

An equation such as equation (A2) occurs in examining airplane stability for the 
controls -fixed case of reference 5 which is modeled in reference 11. Here, r accounts 
for the fact that a vertical gust affects the horizontal tail later than it does the wing. It 
can be shown also that the characteristic quasi -polynomial for the controls -free case of 
reference 5 has the form of equation (4) with 0i = r and 02 = 2r. 

Example 3 : A characteristic quasi -polynomial which occurs in epidemics (ref. 18) is 
given by equation (4) with 

P(s) = T2 (s + a) ^1 = ■’’1 

Q(s) = -1 02 = Ti + T2 

R(s) = 1 
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ROOTS OF QUASI -POLYNOMIAL AS CONTINUOUS FUNCTIONS OF DELAYS 

Theorem : The roots of the quasi -polynomial L(s,0j), where j = 1, 2, . . N, are 
continuous functions of 0j. 

Proof: The quasi -polynomial L(s,0j) is an analytic function, so that its roots are iso- 
lated. Let s* be any one of the roots of L(s,0j) and let p be a positive number 
such that the only root of L(s,0j) contained inside the contour |s - s*| = p is s*. A 
Taylor series expansion of the exponential term e can be used to show that for 
every e > 0, there exists a 6(e) > 0 such that 

L(s,0j) - L(s,0j + A0j)| < e 

for |s - s*i ^ p whenever A0j < 6(e). Choose e< min|L(s,0j)| on js - s*j = p then 
L(s,0j) - L(s,0j + A0j)| < e < |L(s,0j)| 

for A0j sufficiently small and |s - s*| = p. By Rouche's theorem (ref. 19), L(s,0j) 
and L(s,6j + A0j) have the same number of roots inside the circle |s - s*| = p. Since 
p can be made arbitrarily small, it follows that the roots of L(s,0j) must move con- 
tinuously with 0j. 
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partitioning of the delay space by curves generated in the delay space for ^ = 0.2 and various values 
of (Kj,K 2 ) and (Vj,V 2 ). Arrows indicate direction of increasing o). 




)f real part of root -locus curve with respect to for constant 6 
ves of figure 5. Arrows indicate direction of increasing co. 




